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Abstract 

This paper presents an implementation of the approximate statistical moment method for uncertainty propagation 
and robust optimization for a quasi l-D Euler CFD code. Given uncertainties in statistically independent, random, 
normally distributed input variables, a first- and second-order statistical moment matching procedure is performed to 
approximate the uncertainty in the CFD output. Efficient calculation of both first- and second-order sensitivity 
derivatives is required. In order to assess the validity of the approximations, the moments are compared with 
statistical moments generated through Monte Carlo simulations. The uncertainties in the CFD input variables arc 
also incorporated into a robust optimization procedure. For this optimization, statistical moments involving first- 
order sensitivity derivatives appear in the objective function and system constraints. Second-order sensitivity 
derivatives are used in a gradient-based search to successfully execute a robust optimization. The approximate 
methods used throughout the analyses are found to be valid when considering robustness about input parameter 
mean values. 


Nomenclature 

A nozzle area 

a geometric shape parameter 

b geometric shape parameter 

b vector of independent input variables 

F vector of CFD output functions 

g vector of conventional optimization constraints 

k number of standard deviations 

M Mach number at nozzle inlet 

M vector of Mach number at each grid point 
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Minf free-stream Mach number 

Mt target inlet Mach number 

N sample size 

Pb normalized nozzle static back (outlet) pressure 
Q vector of flow-field variables (state variables) 
q mass flux through nozzle 

qt target mass flux through nozzle 

R vector of state equation residuals 

V nozzle volume 

Vt target nozzle volume used for optimization 

x normalized axial position within nozzle 

g standard deviation 

G 2 variance 

superscript: 

mean value 


1 

American Institute of Aeronautics and Astronautics 


Introduction 

Gradient-based optimization of complex 
aerodynamic configurations and their components, 
utilizing high-fidelity Computational Fluid Dynamics 
(CFD) tools, continues as a very active area of research 
(see, for example, Ref. 1 , 2, and 3). In most of the 
CFD-based aerodynamic optimization and design 
studies to date, the input data and parameters have been 
assumed precisely known; we refer to this as 
deterministic or conventional optimization. When 
statistical uncertainties exist in the input data or 
parameters, however, these uncertainties affect the 
design and therefore must be accounted for in the 
optimization. Such optimizations under uncertainty 
have been studied and used in structural design 
disciplines (see, for example, Ref. 4, 5, 6, 7 and 8); we 
refer to these as non-deterministic or robust design 
optimization procedures. 

Sensitivity derivatives (SD) of CFD code output, 
with respect to code input and parameters, contain 
information which can be used to direct the 
optimization search; that is, the objective and constraint 
gradients are functions of the CFD SD. Such SD can 
also be used to accurately approximate the CFD output 
in a small region, such as that near the mean value of a 
random variable. In Ref. 9, it is shown that a statistical 
First Order Second Moment (FOSM) method and 
Automatic Differentiation (AD) can be used to 
efficiently propagate input uncertainties through finite 
element analyses to approximate output uncertainty. 
This uncertainty propagation method is demonstrated 
herein for CFD code. 

An integrated strategy for mitigating the effect of 
uncertainty in simulation-based design is presented in 
Ref. 10; this strategy consists of uncertainty 
quantification, uncertainty propagation, and robust 
design tasks or modules. Two approaches are discussed 
there for propagating uncertainty through sequential 
analysis codes: an extreme condition approach and a 
statistical approach. Both approaches can be efficiently 
implemented using SD. For CFD code, the former 
approach is demonstrated in Ref. 1 1, whereas the latter 
approach is demonstrated herein using second moment 
approximations and SD. These uncertainty propagation 
methods have been developed and are being 
investigated as an alternative to propagation by direct 
Monte Carlo simulation for potentially expensive CFD 
analyses. 

The present paper shows how the approximate 
statistical second moment methods, FOSM and the 
Second Order Second Moment (SOSM) counterpart, 
can be used in conjunction with SD to propagate input 


data uncertainties through CFD code to estimate output 
uncertainties. The FOSM approximation is then used to 
perform sample robust optimizations. For 
demonstration purposes, we assume that the input 
uncertainty quantification is given by independent 
normally-distributed random variables, and we 
demonstrate the strategy of Ref. 10 as applied to a CFD 
code module. This strategy is also applicable to 
correlated and/or non-normally distributed variables; 
however, the analysis and resulting equations become 
more complex. 

The gradient-based robust optimization 
demonstrated herein requires second-order SD from the 
CFD code. A companion paper, Ref. 1 2, presents, 
discusses, and demonstrates the efficient calculation of 
second-order SD from CFD code using a method 
proposed, but not demonstrated, in Ref 13. This 
method, used herein, incorporates FO SD obtained by 
both forward-mode and reverse-mode differentiation in 
a non-iterative scheme to obtain SO SD. 

To date, the only other demonstration or application 
of gradient-based, robust optimization involving 
advanced or high-fidelity (nonlinear) CFD code that we 
have found was just recently presented in Ref. 14 and 
15. The analytical statistical approximation of their 
objective function for robust optimization also required 
second-order SD. However, these studies employed a 
direct numerical random sampling technique to 
compute expected values at each optimization step in 
order to avoid the second-order SD. An example of 
linear aerodynamics involved in multidisciplinary 
performance optimization subject to uncertainty is 
found in Ref. 16. 

Two other aspects need to be pointed out in regard 
to the robust optimization demonstrations for CFD code 
modules presented herein and also in Ref. 14 and 15. 
First, the sources of uncertainty considered were only 
those due to code input parameters involving geometry 
and/or flow conditions; i.e., due to sources external to 
the CFD code simulation. Other computational 
simulation uncertainties, such as those due to physical, 
mathematical and numerical modeling approximations 
(see Ref. 17 and 18) - essentially internal model error 
and uncertainty sources, were not considered. That is, 
the discrete CFD code analysis results were taken to be 
deterministically "certain" herein. Ultimately, all of 
these modeling sources of error and uncertainty must be 
assessed and considered. Sensitivity derivatives can 
also aid in this assessment (Ref. 19) since the adequacy 
of an internal model’s (i.c., algorithm, turbulence, etc.) 
prediction capability generally depends, to some extent, 
on the modeling parameter values specified as input. 
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Second, as discussed in Ref. 15, uncertainty 
classification with respect to an event’s impact (from 
performance loss to catastrophic) and frequency (from 
everyday fluctuation to extremely rare) sets the problem 
formulation and solution procedure. Structural 
reliability techniques typically deal with risk 
assessment of infrequent but catastrophic failure modes, 
identifying the most probable point (MPP) of failure 
and its safety index. Here, we are addressing the 
assessment of everyday operational fluctuations on 
performance loss, not catastrophe. Consequently, we 
are most concerned with aero performance behavior 
due to probable fluctuations, i.e., near the mean of 
probability density functions (pdf). Structural 
reliability assessment is most concerned with 
improbable catastrophic events, i.e., probability in the 
tails of the pdf. Simultaneous consideration of both 
types of uncertainty is discussed in Ref. 16. 

In Ref. 10 an integrated methodology for dealing 
with uncertainty in a simulation-based design is 
proposed and demonstrated for a linkage mechanism 
design. The integrated strategy of Ref. 10 for mitigating 
the effect of uncertainty includes (a) uncertainty 
quantification, (b) uncertainty propagation, and (c) 
robust design. The present study utilizes the strategy 
proposed in Ref 10, but differs in regard to uncertainty 
propagation and application. Here, we are considering 
the influence of uncertainty in CFD code input; that is, 
the effect of uncertainty in input geometry on 
aerodynamic shape-design optimization and the effect 
of uncertainty in flow' conditions on design for flow 
control. 


(FOSM arid SOSM) where the required SD are obtained 
by hand or by AD (see Ref. 9 and 20). The first step in 
both FO and SO analyses is to approximate the CFD 
system output solutions of interest in Taylor series 
form. These approximations are formed to estimate the 
output value for small deviations of the input. 

Given input random variables b={b|,...,b n } with 
mean b = { bj b n } and standard deviations, 
o b = {a b| ,...,a bn ) , the CFD output function, F, firsl- 
and second-order Taylor series approximations are 
given by 


FO: 


F(b)=F(b) + yJ^(b, -b;) 

SO: 

F(b)-F(b)+y|f(b i -b i )+ 


i n n 

+ iy y J_L 
2 't! SfiM-i 


(bj -bj Kbj - bj). 


0 ) 


( 2 ) 


where both first and second derivatives are evaluated at 
the mean values, b . 

One then obtains expected values for the mean (first 
moment) and variance (second moment) of the output 
function, F, which depend on the SD and input 
variances, a b 2 . (Recall the variance is equivalent to the 
square of the standard deviation.) The mean of the 
output function, F, and standard deviation o F , are 
approximated as 


Integrated Statistical Approach 

Our implementation of the three aspects of the 
integrated strategy of Ref. 10 are as follows: 

Uncertainty Quantification 

In this study, we consider the influence of 
uncertainty in CFD input parameterization variables. 
We have assumed that these input variables are 
statistically independent, random, and normally 
distributed about a mean value. This assumption not 
only simplifies the resulting algebra and equations, but 
also serves to quantify input uncertainties. 
Furthermore, it is not an unreasonable assumption for 
input geometric variables subject to random 
manufacturing errors nor for input flow conditions 
subject to random fluctuations. 

Uncertainty Propagation 

Uncertainty propagation is accomplished by 
approximate statistical second moment methods 


FO: 


F = F(b) 



where both first and second derivatives arc evaluated at 
the mean values, b . Note in Eq. (4) that the second- 
order mean output, F , is not at the mean values of 
input b , i.e., F ^ F(b) . 
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Equations (3) represent a FO method and Eq. (4) a 
SO method for examining uncertainty propagation. The 
methods arc straightforward with the difficulty largely 
lying in computation of the SD. The very efficient and 
effective method used here to obtain such derivatives is 
presented in a companion paper, Ref. 12. 

Robust Design 

Conventional optimization for an objective function, 
Obj, that is a function of the CFD output, F, state 
variables, Q, and input variables, b, is expressed in Eq. 
(5). Herein, the CFD state equation residuals, R, are 
represented as an equality constraint, and other system 
constraints, g, are represented as inequality constraints. 
The input variables, b, are precisely known, and all 
functions of b are therefore deterministic. 

min Obj, Obj = Obj(F.Q.b) 

subject to 

R(Q,b) = 0 (5) 

g(F,Q,b) < 0 

For robust design, the conventional optimization, 

Eq. (5), must be treated in a probabilistic manner. 

Given uncertainty in the input variables, b, all functions 
in Eq. (5) are no longer deterministic. The design 

variables are now the mean values, b = { bj b n } , 

where all elements of b arc assumed statistically 
independent and normally distributed with standard 
deviations a b . The state equation residual equality 
constraint, R, is deemed to be satisfied at the expected 
values of Q and b, that is the mean values Q and b for 
the FO approximation. The objective function is cast in 
terms of expected values and becomes a function of F 
and g f . The other constraints are cast into a 
probabilistic statement: the probability that the 
constraints are satisfied is greater than or equal to a 
desired or specified probability, P k . This probability 
statement is transformed (see Ref. 10) into a constraint 
involving mean values and standard deviations under 
the assumption that variables involved are normally 
distributed. The robust optimization can be expressed 
as 

min Obj, Obj = Obj( F ,o Fi Q , b ) 

subject to __ 

R (QJO = 0 (6) 

g(F,Q,b) + ka g <0, 

where k is the number of standard deviations, o g , that 
the constraint g must be displaced in order to achieve 
the desired or specified probability, P k . For the FOSM 
approximation, standard deviations 0 F and a g are of the 
form given in Eq. (3) involving first-order SD. 


Therefore, a gradient-based optimization will then 
require second-order SD to compute the objective and 
constraint gradients. Note that for the SOSM 
approximation, third-order SD would be required for 
these gradients. 

The calculation of second-order SD for CFD code, 
such as those required for SOSM and robust 
optimization with FOSM, was demonstrated in Ref. 13; 
the efficient calculation used herein is demonstrated 
and discussed in a companion paper, Ref. 12. Both 
hand differentiation and AD via the ADIFOR tool (Ref. 
21, 22 and 23) were used. Both conventional and 
robust optimizations were performed using the 
Sequential Quadratic Programming (SQP) method 
option in the Design Optimization Tools, DOT (Ref. 
24). 

Application to Quasi 1-D Euler CFD 

A very simple example has been chosen to 
demonstrate the propagation of input uncertainty 
through CFD code and its effect on optimization. Two 
separate applications are presented; the first involving 
propagation of geometric uncertainties, the second 
involving propagation of flow parameter uncertainties. 
Both uncertainty analyses are performed with quasi 
one-dimensional Euler equations and boundary 
conditions describing subsonic flow through a variable 
area nozzle. The nozzle inlet is located at x = 0 with 
area A(x = 0) = 1 ; the nozzle outlet is at x = 1 . The area 
distribution is given by 

A(x) = 1 - ax + bx 2 . 

The volume, V, occupied by the nozzle, is the 
integration of A(x) over the length x = 0 to x = I 


where a and b are the input geometric parameters. 

Three flow parameters are specified as input 
boundary conditions: the stagnation enthalpy, inlet 
entropy, and outlet static (back) pressure. The quasi 
i-D Euler equation set is symbolically written as the 
state equation in Eq. (5); its residual, R is driven to 
(machine) zero for a solution. 

For supersonic flow through a variable area nozzle, 
shock waves generally appear and the flow solution 
(objective, constraint, etc.) becomes noisy or non- 
smooth (see Ref. 25 and the references cited therein). 
Care must be exercised with respect to obtaining and 
using the SD needed for gradient-based optimization 
(Ref. 25 and 26). Therefore, we chose to bypass issues 
related to this supersonic flow non-smoothness in these 
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initial demonstrations of the present approach for 
uncertainty propagation and robust design for CFD 
code modules. 

Geometric Uncertainty Propagation 

For the discussion of geometric uncertainty 
propagation, geometric shape parameters a and b will 
represent the statistically independent random input 
variables, b. The Mach number distribution through the 
nozzle, M, is viewed here as a component of the state 
variable, Q; its value at the inlet, M, is the CFD output, 
F. Applying the approach previously outlined (recall 
Eq. (3) and (4)) yields the following first- and second- 
order approximations of the output function, M. 


Input random variables: b={a,b) 


CFD output function: F=(M) 


SO: 


M = M(a,b) + 0.5 


k 2 vk 


3a 2 


Om — 


f3M f 


a.. 

V* ", 


3M 

3b 1 


2 ^ - 

Vm 

1 2 

+ 0.5 

3b 2 1 

k 


v J 

f 

3 2 M 

\ 

2 

+ 0.5 

l*' 

/ 


+ ( 10 ) 


+ 0.5 


f d 2 M 2 


3b 


r u b 


Y ^ 2 m 


3a3b 


-a,o 




Predictions of M(a,b), M , and c M for FO (Eq. (7) 
and (9)) and SO (Eq.(8) and (10)) are compared with 
CFD solutions and Monte Carlo analyses based on CFD 
solutions, as given and discussed in the results section. 


Robust Shape Optimiza tion 

Applying the conventional optimization previously 
described yields 


FO Taylor series: min Obj, 

M(a,b) = M(a,b)+ — (a — a) +— (b-b) (7) subject to 

3a 3b 


Obj = Obj(M,a,b) 

R(M,a,b) = 0 
V(a,b) < 0, 


(ID 


SO Taylor series: 


M(a, b) = M(a, b) + — (a - a) + — (b - b) + 
3a 3b 


3 - A/T 

+ — — (a - a)(b - b) + 0.5| 
3a3b 


/ 3 2 M 

3a 2 


(a - a) 2 


\ 


( 8 ) 


+ 0.5 


^ 3“M 
3b 


— (b - b) 


The mean, M , and standard deviation Gm of the 
output function are expressed as 


where the system constraint, V, is a constraint on the 
nozzle volume and depends only on a and b; and our 
objective does not explicitly depend on M. 

Applying the robust optimization previously 
described yields 

min Obj, Obj = Obj( M ,g m , a , b ) 

subject to 

R( M , a , b ) = 0 (12) 

V( a , b )+ ko v ^ 0, 

where 


? I 

3V 

j 

3V V 

o v - = 


J + l 

[db J 


FO: 


M = M(a, b) 


— 


3M 

"sT 




3M 
3b 1 


With a and b subject to statistical uncertainties 
(which may be due to measurement, manufacturing, 
etc.), V becomes uncertain. Since V is linearly 
dependent on a and b, it is also normally distributed. 
Therefore, its standard deviation, Gy, is given exactly 
by Eq. (13). 

To demonstrate the optimizations, a simple target- 
matching problem is selected; a unique answer is 
obtained when an equality volume constraint is 
enforced. The CFD code is run for given a and b; the 
resulting M(a,b) and corresponding V(a,b) are taken as 
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Obj = Obj(M,Minf,Pb) 


the target values Mt and Vt, respectively. For this 
conventional optimization, the objective function and 
constraint function for V of Eq. (11) become 

Obj(M,a,b) = [M(a,b) - Ml] 2 
V(a,b) - Vt = 0 

enforced as ( 1 4) 

V(a,b) - Vt < 0 and Vt - V(a,b) < 0 
for the convenience of the optimizer. 

For robust optimization using the FOSM 
approximation, the corresponding objective and 
constraints on V of Eq. ( 1 2) become 


min Obj 
subject to 

R(M,Minf,Pb) = 0 (16) 

q(Minf,Pb) < 0, 

where q is a constraint on the mass flux through the 
nozzle. 

The robust optimization is expressed as 
min Obj, Obj = Obj( M ,o M , Minf, Pb ) 

subject to 

R(M ,Minf,Pb ) = 0 (17) 

q( Minf, Pb )+ kc q < 0. 


Obj ( M *Gm, a , b ) — [ M (a , b ) - Mt]~ 

V(a , b )-Vt + ka v =0 

similarly enforced as (15) 

V( a , b ) - Vt + ka v < 0 
and Vt- V(a,b)-kc v <0. 

Note that for c a = a b = 0 in Eq. ( 1 5), the 
conventional optimization is obtained. Also, in the 
probabilistic statement of the constraint on V, it is 
assumed that the desired volume is less than or equal to 
Vt. 

Flow Parameter Uncertainty Propagation 

A second example of uncertainty in CFD involves 
fluctuations in input flow parameters. For the 
discussion of flow parameter uncertainty propagation, 
the free-stream Mach number, Minf, and the nozzle 
static back pressure, Pb, will be taken as statistically 
independent random variables. Specifying the free- 
stream Mach number sets the stagnation enthalpy. The 
Mach number distribution through the nozzle, M, is 
again viewed as a component of the state variable, Q; 
its value at the inlet, M, is the CFD output, F. 

Applying the approach previously outlined yields 
equations which are similar to Eq. (7) through (10) but 
with 


For the free stream Mach number, Minf, and the 
nozzle back pressure, Pb, subject to statistical 
uncertainties, the mass flux, q, becomes uncertain. 
Since q is dependent on Minf and Pb, its standard 
deviation, G q , may be approximated by 


r 3q f f 9q ’f 

v ai^F aMinf J + [apb 0ph J ' 


(18) 


Since q is not a linear function of Minf and Pb, Eq. (18) 
is not exact (unlike the previous example where g v was 
exactly known). 

To demonstrate the optimizations, a simple target- 
matching problem is again chosen. The CFD code is 
run for given Minf and Pb; the resulting M and 
corresponding q are taken as the target values Mt and 
qt, respectively. For this conventional optimization, the 
objective function and constraint functions of Eq. (16) 
are 


Obj(M,Minf,Pb) = [M(Minf,Pb) - Mt] 2 
q(Minf,Pb) - qt = 0 

enforced as (19) 

q(Minf,Pb) - qt < 0 and qt - q(Minf,Pb) < 0. 


For robust optimization using the FOSM 
approximation, the corresponding objective and 
constraint on q of Eq. (17) can be shown as 


Input random variables: b={Minf,Pb) 

CFD output function: F={M). 

Again, predictions of M, M , and a M for FO and SO 
approximations are compared with CFD solutions and 
Monte Carlo analyses based on CFD solutions, as given 
and discussed in the next section. 


Obj = Obj( M ,g m , Minf, Pb ) 

= [M(Minf,Pb)-Mt] 2 + o v1 2 
q( Minf, Pb )-qt +kc q =0 

enforced as (20) 

q( Minf, Pb )-qt +ko q < 0 
and qt-q( Minf, Pb )-kc q < 0. 


Robust Design for Flow Control 

The conventional optimization is expressed as 


Again note that for G Min f = G Pb = 0 in Eq. (20), the 
conventional optimization is obtained. Also, in the 
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probabilistic statement of the constraint on q, it is 
assumed that the desired mass flux is less than or equal 
to qt. 


Sample Results & Discussion 

Presentation and discussion of results for the sample 
quasi 1-D Euler CFD problems are divided into four 
topics: function approximations, uncertainty 
propagation, pdf approximations, and robust 
optimization. For the first three topics, the 
approximations are assessed by comparison with direct 
CFD simulations. 

Function Approximations 

It is important to assess the Taylor series output 
function approximations with direct nonlinear CFD 
code simulations prior to presenting uncertainty 
propagation. If the CFD output function, M, is quasi- 
linear with respect to the input variables of interest, one 
can expect First-order approximations to be reasonably 
good; that is, the FO moments given by Eq. (3) should 
match well with the moments produced by a Monte 
Carlo simulation. For a more nonlinear system, one 
naturally expects better accuracy with second order 
approximations; that is, uncertainty analyses which 
include SO terms should yield results which better 
predict the statistical moments produced by the Monte 
Carlo simulation. 

Figures 1-4 show that for F=M(a,b), M behaves as a 
quasi-linear function in the neighborhood of ( a, b ), 
whereas for F=M(Minf, Pb), M is more nonlinear in the 
neighborhood of ( Minf, Pb ). In these Figures, 
approximations of the CFD output functions, M(a,b) 
and M(Minf, Pb), using the first- (FO) and second-order 
(SO) Taylor series (as given in Eq. (7) and (8) for 
M(a,b)), are compared to direct solution of the Euler 
CFD. In each example, two traces were considered 
through the design space. Trace 1 varied the First input 
variable, while the second remained fixed at its mean 
value, and vice versa for trace 2. The required first- 
and second-order SD needed for construction of the FO 
and SO approximations were obtained by hand 
differentiation and AD as discussed and presented in 
Ref. 12. 


nonlinear. At larger deviations from the mean, a linear 
approximation for M(Minf, Pb) loses accuracy. 



Fig. I. Comparison of Function Approximations vs. 
CFD Solution, Input Variable b = b . 



b-b 


Nonlinear behavior of the CFD result is reasonably Fjg 2 Comparison 0 f Function Approximations vs. 

well approximated by the SO result in all plots; CFD Solution, Input Variable a = a . 

however, there does appear to be an inflection point in 
the CFD results given in Fig. 3. Note that the linear FO 
result is a good approximation in the geometric 
example; the flow parameter example is more 
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Minf - Minf 


Fig. 3. Comparison of Function Approximations vs. 
CFD Solution, Input Variable Pb = Pb . 



- 0.30 “- 0.15 0.00 0.15 0.30 

Pb-Pb 

Fig. 4. Comparison of Function Approximations vs. 
CFD Solution, Input Variable Minf = Minf . 

Uncertainty Propagation 

Approximation of the statistical First and second 
moments is done using equations Eq. (9) and (10) 
(geometric example), and corresponding equations for 
the flow parameter example. Again, both first- and 
second-order SD are required and the prediction is 
straightforward, given these derivatives. An 
independent verification of these approximate mean and 
standard deviation values is obtained here using direct 
Monte Carlo simulation with the quasi 1-D Euler CFD 
code and standard statistical analyses of these Monte 


Carlo results. The standard statistical analyses used 
were from MicroSoft ® Excel 2000 and the random 
number generator MZRAN used was from Ref. 27. 
Tables 1 and 2 give results for the mean (first moment) 
and Tables 3 and 4 give results for the standard 
deviation (second moment) value comparisons. The 
input deviations (a a and (7b) or(GM, nt and a Pt> ), are taken 
to be equal and given in the second column of each 
table. The third column in each table gives the result 
from the Monte Carlo simulation, where the sample 
size (N) used was 3,000. The Monte Carlo error in its 
predicted mean is , which is given in the fourth 

column of Tables I and 2, The FO and SO approximate 
predictions are given in the last two columns of each 
table as percent difference from the Monte Carlo 
results. 


Table 1 . Percent Difference from Monte Carlo (MC) 
for FQ and SO Predictions of M (a, b) 



Input o 

M 

% Error 

% difT w/ MC 

% difT w / MC 

Case 

Cj=0|, 

MC 

MC 

FO Predict 

SO Predict 

1 

0.01 

0.4041 

0.0187 

-0.0105 

0.0656 

2 

0.02 

0.4040 

0.0379 

0.0716 

0.1531 

3 

0.04 

0.4054 

0.0756 

-0.2867 

0.0383 

4 

0.06 

0.4055 

0.1142 

-0.3012 

0.4301 

5 

0.08 

0.4096 

0.1557 

-1.3078 

-0.0209 


Table 2. Percent Difference from MC for FO and SO 
Predictions of M (Minf, Pb) 


Case 

Input a 

M 

MC 

% Error 
MC 

%diffw/ MC 
FO Predict 

%diffw/ MC 
SO Predict 

1 

0.0 1 

0.3933 

0.0056 

0.0037 

-0.0269 

2 

0.02 

0.3932 

0.0114 

0.0187 

-0.1034 

3 

0.04 

0.3898 

0.0229 

0.8917 

0.3991 

4 

0.06 

0.3889 

0.0364 

1.1251 

0.0141 


Table 3. Percent Difference from MC for FO and SO 
Predictions of g m . Geometric Example 



Input a 

a M 

%difTw/MC 

%diffw/MC 

Case 

0’a = 0h 

MC 

FO Predict 

SO Predict 

1 

0.01 

0.0102 

-0.5773 

-0.5708 

2 

0.02 

0.0207 

-1.7026 

-1.6769 

3 

0.04 

0.0414 

-1.5794 

-1.4766 

4 

0.06 

0.0625 

-2.2590 

-2.0296 

5 

0.08 

0.0853 

-4.3987 

-4.0001 
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Table 4. Percent Difference form MC for FO and SO 
Predictions of Ov 1; Flow Parameter Example 



Input O 

0\1 

%diffw/ MC 

%diffw/MC 

Case 

OM,nf=C|*i 1 

MC 

FO Predict 

SO Predict 

I 

0.01 

0.0030 

1.1815 

1.2473 

2 

0.02 

0.0062 

-1.5093 

-1.2533 

3 

0.04 

0.0125 

-4.1604 

-3.1680 

4 

0.06 

0.0199 

-4.4070 

-2.1938 


Note that mean and standard deviation 
approximations for output M agree well with the Monte 
Carlo results for small standard deviations in input 
variables, (a a ,a h ) or (a Min f,G Pb ), for both the FO and SO 
predictions. At higher standard deviations, the SO 
results tend to agree better with the Monte Carlo, 
especially in the flow parameter example. Table 3 
illustrates the expected behavior, i.e., SO 
approximations are always better than FO and an 
increase in the standard deviation of input parameters 
monotonically increases the error associated with both 
predictions. Also note that the Monte Carlo simulation 
with a sample size of 3000 limits accuracy of the Monte 
Carlo results. This is apparent in the raggedness of the 
pdfs shown in the following section. 


that in aerodynamic performance optimization using 
CFD, where robustness about the mean is desired, these 
approximations may be good enough. 



Fig. 5. Probability Density Function for M(a,b) for 
a a =a b =0.08. 


(input stdvs=0.02 output std v=0.01 2276) 


Probability Density Function Approximations 

Given a mean and standard deviation of the CFD 
output function (from either a Monte Carlo simulation 
or an FO or SO prediction) and assuming a normal 
distribution, one may then construct a pdf to 
approximate the behavior of the non-deterministic 
output function. This approximation is compared to the 
pdf histogram generated from a Monte Carlo simulation 
in Figs. 5, 6, and 7. The bars depict the Monte Carlo 
histogram, and the solid curve is a normal distribution 
at the Monte Carlo mean value and Monte Carlo 
standard deviation as given in the previous tables. The 
Monte Carlo simulation size of 3000 is certainly not 
sufficient to obtain a smooth pdf. We note that both the 
FO and SO normal distributions are indistinguishable 
from this normal Monte Carlo curve at this scale. It is 
apparent that for either the quasi-linear functional 
dependence on a and b (Fig. 5), or for small input 
standard deviations of the flow parameters (Fig. 6), the 
statistical approximations are good for a significant 
region about the mean but tend to break down in 
predicting the tails of the distribution. This is 
significant, for if one is primarily interested in reliable 
failure predictions, as for structural design, this 
prediction may not be good enough. It is felt, however, 



Fig. 6. Probability Density Function for M(Minf,Pb) 
for 0 Mj n f=Gp b =0.02. 

It is not surprising that a nonlinear CFD output 
function behaves differently for randomness in different 
input variables. For M(Minf, Pb) at higher input 
standard deviations (a M inf=apb =0.06), the pdf of the 
output function is no longer normal. In Fig. 7 one can 
see the non-normal behavior of CFD output given 
normally distributed input variables Minf and Pb. 
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(input stdvs=0.06 output stdv=0.03795) 

0.06 r — — 



Mach Inlet 


Fig. 7. Probability Density Function for M(Minf, Pb) 
for CTM„, t =a H , =0.06. 

Robust Optimizations 

Optimization results were generated using the quasi 
1-D Euler CFD code and the procedure given by Eq. (5) 
and (6). As noted earlier, conventional optimization is 
obtained for a a = G b = 0 or for Omm = G Pb = 0. For the 
FOSM approximation, first-order SD are required to 
obtain a^and G v ando q ; therefore, second-order SD 
will be required for their derivatives, which are used in 
the optimization. The second-order SD were obtained 
by the manner presented in Ref. 12, except for those 
associated with Eq. (18) where the first-order SD were 
finite differenced simply for convenience. 

It is seen from Eq. (1 2) and ( 1 7) that the robust 
optimization results should depend on the probabilistic 
parameters (a a , G b ) or (GMj n} \G Pb ), and k. The desired 
probability, P k , is that from the normal cumulative 
distribution function since g v and G q here are assumed 
to be normally distributed. 

For each robust optimization example, two cases are 
presented. For case 1, P k is fixed at k=l, i.e., 
Pj=84.13%, and the effect of increasing the input 
variable standard deviations is addressed. For case 2, 
the standard deviations of the input variables are fixed 
at 0.01 and P k increases. 

Rob u st Shape Optimization Results - In Table 5, 
results for case 1 of the robust shape optimization are 
displayed. For G a =G b ranging from 0 to 0.08, optimal 
values for the input variables (a, b) are listed. As 
a a =a b increases, so does C\. Accordingly, the mean 
values, (a, b), which minimize the objective function 
and satisfy the probabilistic constraint, become 
increasingly displaced from the target volume, Vt. This 


is shown in Fig. 8^ Mean values (a, b) change, keeping 
the mean value, M(a,b), of the probabilistic output 
near the target value, Mt. The robust design points 
track the dashed curve for M = Mt with some 
displacement due to the C M 2 term of the objective, Eq. 
(15). V (a, b) is displaced from the solid curve V = Vt 
by ka v , as required by the probabilistic constraint. This 
displacement can be viewed as the solution dependent 
or "effective" safety factor. 


Table 5. Robust Shape Optimization Results with 
Increasing Input Parameter a for k-1 


a ;i =a h 

a 

b 

Obj 

M 

CTm 

Ox 

0.00 

0.6001 

0.3001 

0.0000 

0.4043 

0.0000 

0.0000 








0.04 

0.7338 

0.4286 

0.0016 

0.4018 

0.0406 

B9 

0.06 





0.0607 

ESI 

0.08 

0.8534 

0.5358 

0.0065 

0.3941 

0.0804 

0.0480 


In Fig. 9 the changing area distribution of the robust 
optimization is illustrated. As the standard deviations 
of design variables (a,b) increase, the optimal non- 
deterministic volume, V (a, b), significantly decreases. 



Fig. 8. Optimization results in design space (a,b), 
P k fixed at Pj. 



Fig. 9. Nozzle area distributions, P k fixed at P h 
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The results for case 2 of the robust shape 
optimization, where o a =o b is fixed at 0.01, and P k 
increases from 50 percent to 99.99 percent (k=0 to 4) 
are given in Table 6. Again mean values (a, b) change, 
keeping the mean value, M(a,b), of the probabilistic 
output near the target value, Mt. Since a a -G b remains 
small, the g m “ term of the objective remains small, and 
the displacement of M from the dashed line depicting 
Mt due to the G M ‘ term remains small as shown in 
Fig. 10. With an increase in P k , V (a, b) is displaced 
from the solid curve V = Vt by ka v , as required by the 
probabilistic constraint. Accordingly, the mean values, 
(a,b) , which minimize the objective function and 
satisfy the constraint, again become increasingly 
displaced from those at the target volume, Vt. Note the 
significant displacement of the solution from the target 
volume when P k is large, i.e., when one is attempting to 
incorporate the tails of the pdf. In order to increase the 
probability of constraint satisfaction from 97.77 percent 
to 99.99 percent, one sees a significant change in (a, b) 
for a mere gain of 2 percent in constraint satisfaction. 


Table 6. Robust Shape Optimization Results with 
Increasing P k for q=0.01 


K 


a 

b 

Obj 

M 


Ov 

0 

0.5000 

0.5996 

0.2995 

0.000104 

0.4041 

0.0101 

0.006 

1 

0.8413 

0.6246 

0.3189 

0.000118 

0.4004 

_0.010l_ 

0.006 

2 

0.9772 

0.6698 

0.3687 

0.000104 

0.4041 

J).0101 

0.006 

3 

0.9986 

0.7052 

0.4037 

0.000104 

0.4042 

0.0102 

0.006 

4 

0.9999 

0.7406 

0.4388 

0.000104 

0.4043 

0.0102 

0.006 



o fixed at 0.01 . 

Robust Design for Flow Control Results - Similar 
results are seen in the flow parameter example. In 
Table 7, the results for case 1 are displayed. For 


^Mim = CT|> b ranging from 0 to 0.08, optimal values for 
the input variables (Minf,Pb) are listed. As GMinf^ G Pb 
increases, so does G q . Accordingly, the mean values, 
(Minf,Pb), which minimize the objective function and 
satisfy the constraint, become increasingly displaced 
from the target mass flux, qt. This is shown in Fig. 1 1 . 
Mean values (Minf,Pb) change, keeping the mean 
value, M(Minf,Pb), of the probabilistic output near the 
target value, Mt. The robust design points track the 
dashed curve for M = Mt with displacement due to the 
o M 2 term of the objective, Eq. (17). The optimized 
mass flux, q (Minf,Pb), is displaced from the solid 
curve q = qt by ko M as required by the probabilistic 
constraint. 


Table 7. Robust Design for Flow Control Results with 
Increasing Input Parameter g for k=l 


O'Min — Opb 

Mini 

Pb 

Obj 

M 

Om 


0 

0.3000 

0.8000 

0.0000 

0.3933 

0.0000 

0.0000 

0.02 

0.2861 

0.7883 

0.0001 

0.3974 

0.0116 

0.0058 

0.04 

0.2655 

0.7801 

0.0005 

0.3985 

0.0231 

0.0112 

0.06 

0.2555 

0.7653 

0.0012 

0.4050 

0.0327 

lo.OI63_ 

0.08 

0.2468 

0.7498 

0.0020 

0.41 18 

0.0407 

0.0209 



Free- Stream Mach No. (Mlnf) 

Fig. 1 1. Optimization results in design space 
(MinfiPb), P k fixed at Pj. 


The results for case 2 of the robust design for flow 
control, where a Mmt =G Pb is fixed at 0.01 , and P k 
increases from 50 percent to 99.99 percent, (k=0 to 4) 
arc given in Table 8. Again, mean values 
(Minf,Pb) change, keeping the mean value, 
M(Minf,Pb), of the probabilistic output near the target 
value, Mt. As in the preceding example, since 
tfMini-tfpb remains small, the c M ~ term of the objective 
remains small and the displacement due to the Gm" term 
remains small, as shown in Fig. 1 2. With an increase in 
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P k , q (MinfiPb) is displaced from the solid curve q = qt 
by ka q , as required by the probab[listi£ constraint. 
Accordingly, the mean values, (MinfiPb), which 
minimize the objective function and satisfy the 
constraint again become increasingly displaced from 
the target mass flux, qt. Again, note the significant 
displacement from the target mass flux incurred in the 
higher probability optimizations, i.e., when one is 
attempting to incorporate the tails of the pdf. 


Table 8. Robust Design for Flow Control Results with 
Increasing P k for a = 0.01 


k 

Pk 

Mint 

Pb 

Obj 

M 

Gm 

G q 

_0 

1 

0.5000 

0.3000 

0.8000 

0.00003 

0.3933 

0.0060 

0.0030 

0.8413 

0.2919 

0.7953 

0.00003 

0.3945 

0.0059 

0.0029 

2 

0.9772 

0.2825 

0.7916 

0.00003 

0.3949 

0.0059 

0.0029 

3 

0.9986 

0.2688 

JX7896 

0.00003 

0.3936 

0.0060 

0.0028 

4 

0.9999 

0.2598 

0.7867 

0.00003 

0.3938 

0.0060 

0.0028 



Free- Stream Mach No. (Mint) 

Fig. 12. Optimization results in design space 
(Minf,Pb), a fixed at 0.0 1 . 


Concluding Remarks and Challenges 

The present results represent an implementation of 
the approximate statistical moment method for 
uncertainty propagation and robust optimization for a 
quasi 1-D Euler CFD code. Assuming statistically 
independent, random, normally distributed input 
variables, a first- and second-order statistical moment 
matching procedure was performed to approximate the 
uncertainty in the CFD output. Efficient calculation of 
both first- and second-order sensitivity derivatives was 
employed and the validity of the approximations was 
assessed by comparison with statistical moments 
generated through Monte Carlo simulations. The 


uncertainties in the CFD input variables were 
incorporated into a robust optimization procedure 
where statistical moments involving first-order 
sensitivity derivatives appeared in the objective 
function and system constraints. Second-order 
sensitivity derivatives were used in a gradient-based 
robust optimization. The approximate methods used 
throughout the analyses were found to be valid when 
considering robustness about input parameter mean 
values. 

Collectively, these results demonstrate the 
possibility for an approach to treat input parameter 
uncertainty and its propagation in gradient-based design 
optimization that is governed by complex CFD analysis 
solutions. It has been demonstrated on a very simple 
CFD code and problem; there are computational 
resource issues to be addressed in application to 
significant 2-D and 3-D CFD codes and problems. 

Some of these are addressed in the companion paper. 
Ref. 12. Work is presently in progress regarding 
application and demonstration using 2-D Euler CFD 
code. 
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